Thermally driven vertical convection (VC) – the flow in a box heated on one side and cooled on the other side, is investigated using direct numerical simulations with Rayleigh numbers over the wide range of $10^7\le Ra\le 10^{14}$ and a fixed Prandtl number $Pr=10$ in a two-dimensional convection cell with unit aspect ratio. It is found that the dependence of the mean vertical centre temperature gradient $S$ on $Ra$ shows three different regimes: in regime I ($Ra \lesssim 5\times 10^{10}$), $S$ is almost independent of $Ra$; in the newly identified regime II ($5\times 10^{10} \lesssim Ra \lesssim 10^{13}$), $S$ first increases with increasing $Ra$ (regime $\textrm {{II}}_a$), reaches its maximum and then decreases again (regime $\textrm {{II}}_b$); and in regime III ($Ra\gtrsim 10^{13}$), $S$ again becomes only weakly dependent on $Ra$, being slightly smaller than in regime I. The transition from regime I to regime II is related to the onset of unsteady flows arising from the ejection of plumes from the sidewall boundary layers. The maximum of $S$ occurs when these plumes are ejected over approximately half of the area (downstream) of the sidewalls. The onset of regime III is characterized by the appearance of layered structures near the top and bottom horizontal walls. The flow in regime III is characterized by a well-mixed bulk region owing to continuous ejection of plumes over large fractions of the sidewalls, and, as a result of the efficient mixing, the mean temperature gradient in the centre $S$ is smaller than that of regime I. In the three different regimes, significantly different flow organizations are identified: in regime I and regime $\textrm {{II}}_a$, the location of the maximal horizontal velocity is close to the top and bottom walls; however, in regime $\textrm {{II}}_b$ and regime III, banded zonal flow structures develop and the maximal horizontal velocity now is in the bulk region. The different flow organizations in the three regimes are also reflected in the scaling exponents in the effective power law scalings $Nu\sim Ra^\beta$ and $Re\sim Ra^\gamma$. Here, $Nu$ is the Nusselt number and $Re$ is the Reynolds number based on maximal vertical velocity (averaged over vertical direction). In regime I, the fitted scaling exponents ($\beta \approx 0.26$ and $\gamma \approx 0.51$) are in excellent agreement with the theoretical predictions of $\beta =1/4$ and $\gamma =1/2$ for laminar VC (Shishkina, Phys. Rev. E., vol. 93, 2016, 051102). However, in regimes II and III, $\beta$ increases to a value close to 1/3 and $\gamma$ decreases to a value close to 4/9. The stronger $Ra$ dependence of $Nu$ is related to the ejection of plumes and the larger local heat flux at the walls. The mean kinetic dissipation rate also shows different scaling relations with $Ra$ in the different regimes.